WGS分析笔记（3）- bam文件的处理

• 上一篇记录了mapping这一步的软件选择，也讲到了对于sam文件如何考虑MAPQ的过滤，这篇我主要想记录一下bam文件在进行call variation之前的处理。
  

• 包括MAPQ的筛选等都是这个处理的一部分。
  

• 既然要处理bam文件，不得不提bam文件的格式。
  

• 处理生物信息的数据时会发现，文件有各种各样的格式，眼花缭乱，这也是我一开始接触课题接触分析流程时的感受。但是多和这些文件打交道以后，会发现，大多数的不同格式的文件其本质都是文本文件，为了某种特殊的处理或记录需要，按一定的规则记录信息。包括之前接触的fasta、fastq文件，也包括sam文件，以后后面步骤会接触到的vcf文件等。
  

• bam文件是sam文件的二进制格式，这里贴一张图，来说明一下为什么要转换成bam文件来处理。
  
  

• 可以看到，不经过处理的sam文件是非常大的，非常不利于存储和处理，而转换格式后的bam文件小很多，所以很多处理软件也是针对bam文件进行开发的。
  

• 包括bwa，bowtie2的输出文件，都是sam文件，但是sam文件具体是什么样子的，我这里不会展开来讲，一是因为官网说明书已经说的很清楚了，二是因为你随手一个百度、谷歌（如果你翻得动墙）就有很多很多的人发帖介绍，而内容大体都是类似的，我自认为也没有什么能比他们讲得更好的。但是了解这个文件的内容确实是非常重要的。
  

• 以@开头的行记录了header信息，之后的行记录了每个reads的mapping信息，一般对于这样的sam文件我们要做的处理大体就是先转换为bam文件，进行sort，再进行merge，最后mark duplicates，并建立索引。接下来我将一步步介绍怎么去处理，以及为什么我是这么处理的。
  

• 在开始之间先使用samtools对参考序列进行索引建立，作用是用于samtools软件快速识别，这一步也是一劳永逸的，只要不把索引文件删除。
  

  $ samtools faidx re.fa
  # 得到索引文件：re.fa.fai

• 首先是把sam文件转换成bam文件，我们通过samtools view来实现，代码如下：
  

  $ samtools view -S in.sam -b > out.bam
  $ samtools view -b in.bam -S > out.sam

• 但实际上，在真正的操作中，我们是不会保留sam文件的，甚至是不会产生sam文件的，因此，我们通常这么来写命令。
  

  $ bwa mem -t 12 -M -R "@RG\tID:W2018001\tPL:ILLUMINA\tLB:W2018001\tSM:W2018001" /your/path/of/reference/ucsc.hg19.fasta 1.fq.gz 2.fq.gz | samtools view -q 1 -Shb - > W2018001.bam
  # 关于“|”之前的我就不解释了，看我上一篇简书
  # -q 1 ：筛选MAPQ用的，意为保留MAPQ >= 1的记录，上一篇简书中讨论过关于MAPQ的由来和分布，这里用“1”也是保险起见，但实际上没啥区别，后期处理上（vcf处理），我们会以更严格的阈值去过滤
  # -h ：表示保留header信息
  # -S，-b：表示格式，S指的是sam格式，b指的是bam格式

• 这样操作的好处是直接可以得到bam文件，不必占用大量的空间存储sam文件。如果你实在是需要sam文件也可以转换回去，但如果你只是想查看一下，也可以通过以下方式实现。还是很方便的。
  

  $ samtools view -h xxx.bam | less

• 在测序的时候序列是随机打断的，所以reads也是随机测序记录的，进行比对的时候，产生的结果自然也是乱序的，为了后续分析的便利，将bam文件进行排序。事实上，后续很多分析都建立在已经排完序的前提下。关于排序可以通过以下命令完成。
  

  $ samtools sort -@ 6 -m 4G -O bam -o sorted.bam W2018001.bam
  # @：线程数
  # m：每个线程分配的最大内存
  # O：输出文件格式
  # o：输出文件的名字
  # 输入文件放在最后

• 接下来要做的是merge，这个不是所有的样本都需要做的！！！
  

• 之前有提到，WGS的数据比较大，对于一个样本可能有多对文件，一般有两个思路，一个是先对原始的fastq文件进行merge，一个是对后面的bam文件进行merge。那么我选用的是后者。实现脚本如下：
  

  $ samtools merge -@ 6 -c sorted.merged.bam *.sorted.bam
  # @：线程数
  # c：当输入bam文件的header一样时，最终输出一个header
  # 当然也可以用-h file，定义你想要的header，file里的内容就是sam文件header的格式
  # 第一个bam是输出的文件名，后面排列输入的文件名，我这里用了通配符‘*’，要保证该目录下‘.sorted.bam’结尾的都是你要输入的文件
  # 当然也可以用-b file，file文件里罗列要merge的文件名，一行一个文件名

• 下一步就是remove duplicates，为什么要进行这一步呢？先来说一下测序，我们都知道人的基因组是很庞大的（约30亿个碱基对），在测序的时候，先会把基因组的DNA序列通过超声震荡随机打断成150bp的片段，那么从概率上来讲，出现同样的片段（开始和结束位置都一样）的几率是极小的。但是由于PCR对某些位置有偏好的扩增，会导致一些一模一样的reads存在。这些reads其实是一个片段的扩增导致的，多出来的reads，对该区段突变的判断是没什么贡献的，如果不加处理，反而会大大增加那个位置的深度，导致某些结果的不准确。
  

• 如下图所示，箭头所指的reads，就是duplicates，我们一般采取的策略是标记或者去除，这样的话，下一步call variation的时候会不考虑这些reads。
  
  

• 关于这一步，有很多软件可以实现，比较多用的是picard和samtools rmdup。我这里用的是GATK包里集成的picard的MarkDuplicates。关于picard和samtools rmdup的效果其实大家可以自己试一下，我很早之前做过的试验是，samtools rmdup速度很快，但是去除的效果稍差。大家用的最多的也是picard。
  

  $ java -jar /you/path/of/gatk/gatk-package-4.0.10.1-local.jar \
  	MarkDuplicates \
  	-INPUT sorted.merged.bam \
  	-OUTPUT sorted.merged.markdup.bam \
  	-METRICS_FILE markdup_metrics.txt
  # 也可以加上“REMOVE_DUPLICATES=true”来去除掉这些duplicates，不然就只是标记

• 到了这一步基本上就处理的差不多了，可以进行call variation了，但是这里还有一步建立索引，这十分的重要！！！！
  

• 必须对bam文件进行排序后，才能进行index，否则会报错。建立索引后将产生后缀为.bai的文件，用于快速的随机处理。很多情况下需要有bai文件的存在，特别是显示序列比对情况下。建立索引很简单。
  

  $ samtools index sorted.merged.markdup.bam

• 到了这一步就基本上完事了，可以通过IGV或tview来查看情况，这都需要建立索引，且索引文件和bam文件在同一个目录下。
  

  $ samtools tview sorted.merged.markdup.bam re.fa
  # 可以用-p chr:pos（加在tview和sorted.merged.markdup.bam之间）指定要查看的位置
  # 也可以进去后用敲‘g’输入`chr10:10,000,000' 或 `=10,000,000'查看指定的位置，敲'?'了解更多说明，q退出

• 下图就是效果了，用“？”，打开左边的帮助界面，其中圆点表示正链比对，逗号表示负链比对，星号表示插入。不同的颜色代表不同的含义，具体怎么调看帮助框。要是觉得不好看的话可以用IGV查看。IGV的效果就是上图duplicates的效果。
  
  

• 同时对于得到的bam文件也可以进行一下查看，对于bam的统计软件就更多了，我这里网罗了一些帖子上的推荐以及我自己日常使用的软件，感兴趣的可以自己去下载来使用一下。
  

###samtools

• 这是个强大的软件，也自带一些统计工具，上篇简书其实我们就用到了，主要是四个：idxstats，depth，stats，flagstat
  

  $ samtools depth sorted.merged.markdup.bam
  	示例结果
  		#chr    pos depth
  		chr1    1   5
  	结果可以得到染色体名称、位点位置、覆盖深度
  	-a:输出所有位点，包括零深度的位点
  	-a -a --aa:完全输出所有位点，包括未使用到的参考序列
  	-b FILE：计算BED文件中指定位置或区域的深度
  	-f FILE：指定bam文件
  	-l INT:忽略小于此INT值的reads
  	-q INT：只计算碱基质量大于此值的reads
  	-Q INT：只计算比对质量大于此值的reads
  	-r CHR:FROM-END：只计算指定区域的reads
  $ samtools idxstats sorted.merged.markdup.bam  #需要预先进行sort和index
  	示例结果
  		#ref    sequence_length mapped_reads    unmapped_reads
  		chr1    195471971    6112404    0
  	结果可看出，分别统计染色体名称、序列长度、mapped数目，以及unmapped数目
  $ samtools flagstat sorted.merged.markdup.bam
  	示例结果：
  		20607872 + 0 in total (QC-passed reads + QC-failed reads) #总共的reads数
  		0 + 0 duplicates #重复reads的数目
  		19372694 + 0 mapped (94.01%:-nan%) #总体上reads的数目以及匹配率；可能有有小偏差
  		20607872 + 0 paired in sequencing  #paired reads的数目
  		10301155 + 0 read1 #reads1的数目
  		10306717 + 0 read2 #reads2的数目
  		11228982 + 0 properly paired (54.49%:-nan%) #完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
  		18965125 + 0 with itself and mate mapped#两条都比对到参考序列上的reads数，但是并不一定比对到同一条染色体上
  		407569 + 0 singletons (1.98%:-nan%)#只有一条比对到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
  		3059705 + 0 with mate mapped to a different chr#两条分别比对到两条不同的染色体的reads数
  		1712129 + 0 with mate mapped to a different chr (mapQ>=5)#两条分别比对到不同染色体的且比对质量值大于5的数量
  	# 说明：
  		1.bwa的mem比对方法生成的bam文件保留sencondly比对的结果。所以使用flagstat给出的结果会有偏差。
  		2.每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以“PASS + FAILED”格式显示
  $ samtools stats sorted.merged.markdup.bam
   	# 对bam文件进行详细的统计，也可只统计某一染色体的某一区域chr：from-to
  	结果包括：
  		Summary Numbers,raw total sequences,filtered sequences, reads mapped, reads mapped and paired,reads properly paired等信息
  		Fragment Qualitites #根据cycle统计每个位点上的碱基质量分布
  		Coverage distribution #深度为1，2，3，，，的碱基数目
  		ACGT content per cycle #ACGT在每个cycle中的比例
  		Insert sizes #插入长度的统计
  		Read lengths #read的长度分布
  	stats出来的结果可以使用plot-bamstats来画图(samtools目录下misc目录中)
  	$ plot-bamstats -p outdir/ sorted.merged.markdup.bam.stats

• 下图就是plot-bamstats操作的输出结果了，可以看到有很多的图。效果还是很好的。
  
  

• 更多关于samtools的工具以及上文提到的工具的其余参数请参考官网：http://www.htslib.org/doc/samtools.html
  

###RSeQC

• 这也是上篇中提到过的，所以安装方式和使用就不赘述了，其实里面还有一些其余实用的工具，当然这款软件的最初使用对象是RNA-seq，但并不影响使用，有些工具是通用的，有一点要注意的是，bam_stat.py里的unique mapping的默认阈值是MAPQ>=30，当然可以通过参数来修改，这个在结果的理解上大家要注意。
  

###bedtools

• 这是一个经常使用也很实用的软件，我经常会用来截取bam然后在igv上看突变的情况，师姐推荐其中的mutlicov进行覆盖度的统计。我粗略的看了一下bedtools的说明书，用于coverage统计的应该还有coverage，genomecov。感兴趣的可以尝试一下。
  

  bedtools:
  	$ wget https://github.com/arq5x/bedtools2/releases/download/v2.27.1/bedtools-2.27.1.tar.gz
  	$ tar -zxvf bedtools-2.27.1.tar.gz
  	$ cd bedtools2
  	$ make
  # 脚本在bin/下的bedtools
  # Ubuntu用户也可以使用下述命令来下载：
  $ sudo apt-get install bedtools

• 截取bam文件，查看igv可以用以下命令：
  

  $ bedtools intersect -a sorted.merged.markdup.bam -b region.sorted.bed > target_reads.bam && samtools index target_reads.bam 
  #其中bed文件的格式可以参考：
  #染色体号  起始位置  终止位置
  chr1	226173187	226173247
  #用"\t"分隔

###GATK

• GATK不仅可以call variation，里面还包含了很多其他用途的工具包，其中有一个工具叫DepthOfCoverage，可以统计depth和coverage，但是在panel中比较适用，因为有bed范围，且比较小。这个工具的速度是比较慢的，要在全基因组上做不太现实。所以我本人也没去使用。
  

###BAMStats

• 一款比较早的bam统计软件，没用过，但是看说明使用是极其简单了，说一下怎么安装。感兴趣的可以自己试一下，很简单。
  

  $ wget https://jaist.dl.sourceforge.net/project/bamstats/BAMStats-1.25.zip
  $ unzip BAMStats-1.25.zip

###bamdst

• 一款简单好用的深度统计软件。
  

  $ git clone https://github.com/shiquan/bamdst.git
  $ cd bamdst/
  $ make

• 安装好后，需要准备.bed文件及.bam文件，以软件提供的MT-RNR1.bed和test.bam为例，使用命令如下：
  

  $ bamdst -p MT-RNR1.bed -o ./ test.bam

• 输出的结果包含7个文件，为：
  

  • -coverage.report
  • -cumu.plot
  • -insert.plot
  • -chromosome.report
  • -region.tsv.gz
  • -depth.tsv.gz
  • -uncover.bed

• 主要看一下coverage.report文件，里面包含了大量信息。
  

###qualimap

• 这算是压轴了吧，这个软件是我师姐推荐的，安装使用都比较容易，给出的是PDF或HTML的报告，报告中的信息特别丰富，还有一堆的图，所以在我自己的脚本中也会对每个样本使用该软件统计，简述一下安装和使用。
  

  # 第一步：下载
  $ mkdir qualimap
  $ cd qualimap
  $ wget https://bitbucket.org/kokonech/qualimap/downloads/qualimap_v2.2.1.zip
  $ unzip qualimap_v2.2.1.zip
  $ cd qualimap_v2.2.1
  # 第二步安装相关的软件
  # java
  # 这个应该都有，这里就是贴一下官网的步骤
  $ sudo apt-get install openjdk-6-jre
  # R
  # 官网上也有，我贴的是自己以前安装的记录
  $ apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9
  $ apt-get update
  $ apt-get install r-base
  $ apt install r-base-core
  # R包安装，两个方法，第一个听说容易报错，我用的是第二个
  $ Rscript scripts/installDependencies.r
  # 或
  $ R
  > install.packages("optparse")
  > source("http://bioconductor.org/biocLite.R")
  > biocLite(c("NOISeq", "Repitools", "Rsamtools", "GenomicFeatures", "rtracklayer"))
  > q()

• 使用也简单，主要分为带gtf文件和不带gtf文件，全基因组的话一般不带gtf文件，然后bamqc其实也只是这个软件中的一个工具，其他的工具感兴趣的也可以看看。
  

  $ qualimap bamqc -bam sorted.merged.markdup.bam  --java-mem-size=20G -c -nt 16 -outdir bamqc -outfile bamqc.pdf -outformat PDF:HTML
  # 参数也没有什么特别要注意的，基本上默认的就可以

• 找了一个示例结果，发现有23页，我这里就不贴了，大家可以自己尝试一下。
  

• 最后贴两张图，是我自己写的脚本得到的深度分布，累积曲线以及覆盖率，因为是自己写的，所以比较糙，横纵坐标标题什么的一律没写。
  
  

• 上图可以看到，深度分布还是比较正态的，最多的30X左右，至于左边0X为什么这么高，是因为参考基因组有些位置就是N，还有一些位置就是我的样本没有覆盖到。
  

• 上图可以看到，小于10X的数据的不到2%，超过50%的数据都是高于30X的，还是不错的。
  

• 上图我按不同的染色体进行统计覆盖率，去掉了其余的一些未知染色体位置的序列上的信息（这个信息具体要了解参考基因组的特性，比如有些序列目前能明确在人身上，却不知道具体在哪个染色体等，这些信息也是包含在参考基因组中的，仔细看参考基因组会发现，除了22条常染色体，2条性染色体，还有很多其他的序列），这个覆盖率并不是我想想的那般全体高于99%，也没公司给的报告描述的那么好，主要是因为MAPQ的过滤导致了这样的结果，但是总体的覆盖率还可以。

• 水平有限，要是存在什么错误请指出，可发送邮件至shiyuant@outlook.com！请大家多多批评指正，相互交流，共同成长，谢谢！！！